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ABSTRACT 

We present a new Particle-Mesh cosmological N-body code for accurately solving the mod- 
ified Poisson equation of the Quasi Linear formulation of MOND. We generate initial con- 
ditions for the Angus (2009) cosmological model, which is identical to AC DM except that 
the cold dark matter is switched for a single species of thermal sterile neutrinos. We set the 
initial conditions at z — 250 for a (512 Mpc/Kf box with 256'' particles and we evolve them 
down to z = 0. We clearly demonstrate the necessity of MOND for developing the large scale 
structure in a hot dark matter cosmology and contradict the naive expectation that MOND 
cannot form galaxy clusters. We find that the correct order of magnitude of X-ray clusters 
(with Tx > 4.5 keV) can be formed, but that we overpredict the number of very rich clusters 
and seriously underpredict the number of lower mass clusters. The latter is a shortcoming of 
the resolution of our simulations, whereas we suggest that the over production of very rich 
clusters might be prevented by incorporating a MOND acceleration constant that varies with 
redshift and an expansion history that cannot be described by the usual Friedmann models. 
We present evidence that suggests the density profiles of our simulated clusters are compat- 
ible with those of observed X-ray clusters in MOND. It remains to be seen if the low mass 
end of the cluster mass function can be reproduced and if the high densities of dark matter in 
the central 20 kpc of groups and clusters of galaxies, measured in the MOND framework, can 
be achieved. As a last test, we computed the relative velocity between pairs of halos within 
10 Mpc and find that pairs with velocities larger than 3000 kms^^, like the bullet cluster, 
can form without difficulty. 
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1 INTRODUCTION 

Cosmological simulations and other calculations of the for- 
mation of large scale structure in the Modified Newtonian Dynam- 
ics (MOND) o f iMilgronj fl983h have a substantial literature, de- 
spite the absence of an agreed upon cosmological model. Early 
work by ISandersI ( Il998h crucially demonstrated that the represen- 
tation of the expansion of the Universe by a scale factor, derived 
from the Friedmann-Robertson-Walker (FRW) metric, is justified 
in MOND becau se expansion w ill indeed be uniform, contrary to 
earlier claims by iFelteiil ( Il984l) . In addition to this, the results of 
Big Bang Nucleosynthesis (assuming the number of neutrinos is 
not modified) were shown to be unaltered by MOND, and the natu- 
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ral result that structure formation proceeds more rapidly in MOND 
was verified. 

Sanders (2001') continued this theme by modifying the CMB- 
FAST code of Seljak & Zaldaniaga ( 1996) for a two field MOND- 
like Lagrangian based theory and showed that, by and large, the 
observed distribution of galaxies (the power spectrum at the time) 
appeared to be reproduced in a low density, baryon only, vacuum 
energy dominated model. Unfortunately for that model, today's 
precision cosmology contradicts it at high significance, but never- 
theless, a solid basis was set which proved matter over-densities, 
which evolve ponderously in Newton's gravity, can evolve with 
great vigour in Milgrom's modified dynamics (see his Fig 2). 

Following these pivotal works, iNusserl ( |2002|) made the 
first foray into MOND cosmological simulations by modifying 
a particle-mesh (PM) code to include an approximation to the 
MOND equation and then performed a series of cosmological sim- 
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ulations. Unfortunately, the simulations were run without any vac- 
uum energy, but generally showed the alarming trend for MOND 
to produce too much structure on the scales A: = 0.1 — 1.0 k/Mpc. 
This conclusion prevailed even in an open Universe with f2m=0.03 
and the acceleration constant of MOND being reduced to 1/12 
of its typical value. One thing that is uncertain is whether or 
not the initial conditions were self-consistent for a distribution of 
baryons in a low density Universe. For instance, the power spec- 
trum at redshift z = 100 for k > 0.5h/Mpc should be less than 
W^^'^ k'^ /Mpc^ , in contrast to a cold dark matter dominated Uni- 
verse (f2™^^J^0)Jorwhich^P(fc) > 10"^ for all k < i.Qh/Mpc. 

iKnebe & Gib"sor] j2004l) incorporated the algebraic MOND 
equation into their Multi Level Adaptive Particle Mesh code (called 
MLAPM) and, using it, made some similar comparisons to the 
AC DM model with MOND as lNussejjioO^) 

In a seminal piece of work, Skordis et al. ('2006') linearised the 
equations of motion from Bekenstein ( 2004) and incorporated them 
into CMBFAST to derive the linear evolution of the formation of 
structure in that theory. They showed that, with the addition of 3 
species of 2.75 eV active neutrinos (and a substantially increased 
MOND acceleration constant), the galaxy power spectrum, as mea- 
sured by the SDSS survey, could be matched by the evolved power 
spectrum from the code. Whether galaxies or clusters would go on 
and form the correct distr ibution is unknown, b ut highly doubtful, 
because shortly thereafter lDodelson & Liguoril (2006) perfomed a 
similar calculation and showed that in the framework of TeVeS, 
there was no way to form the coiTect distribution of cosmic struc- 
ture without non-baryonic dark matter. 

Leaving aside relativistic counterparts of MO ND, a crucial 
development on the MLAPM co de was made by iLlinares et al.l 
( l2008h and further used in ILlinares et al.l ( |2009|) because, instead 
of "merely" solving the algebraic MOND equation, they solved the 
mo mentum and energy conservin g modified Poisson equation of 
the lBekenstein & Milgrom' ('1984^) theory using the same technique 
as originally outline d by Brada & Milgrom ( 1999) and further used 
bv lXiret & CombesI (|2007|) for galaxy simulations. 

Although incredibly important stepping stones, the obvious 
flaw with all these previous works is that, although they solve some 
variant on the MOND equation in their simulations, they start from 
arbitrary initial conditions. None were derivatives of models that 
were consistent with observations of the anisotropies in the cos- 
mic microwave background (CMB) nor the presence of dark mat- 
ter in c l usters of galaxies. This no longer has to be the case because 
lAngu j (12009') has proposed a model that can allow MOND to sat- 
isfy many cosmological constraints. 

The model currently employs a simple ansatz - that whichever 
relativistic version of MOND is cortect, the FRW metric is still the 
correct description of the expansion history of the Universe, which 
is a perfectly reason able assumption expected from relativistic ana- 
logu es of MOND (IFeix et"ani2010l : Iciifton & Zlosni3l2010l. see 
also ISkordis & Zlosnikl 201 ll) . Given that typical accelerations at 
high z are of the order of the MOND acceleration constant, this is a 
justified simplification at those redshifts. Even were it not, there is 
no reason to believe that the MOND acceleration constant, which is 
suggested to be linked to the presence of dark energy, need be con- 
stant with redshift. There is, however, the very logical possibility 
that the expansion history cannot be described by a standard Fried- 
mann model and this is a topic that is addressed further in i]5.1l 

The idea behind the model is that MOND has been argued 
to be consistent with the dynamics of all galactic and sub-galactic 
dynamics, i ncluding the d warf spheroidal galaxies surrounding the 
Milky Way l lAnguj|200§) and the recently analysed THINGS ro- 



tation curves ( iGentile et al.ll201 ll) . without the need for cold dark 
matter. Although there are still contentious issues surrounding these 
systems, the balance of evidence suggests MOND can at least do 
as good a job of fitting their dynamics as a cold dark matter halo 
which has several extra free parameters. When th is is added to the 
evidence from the baryonic TuUy-Fisher relation (lMcGaughll201 ih 
it provides a strong case that MOND is a good description of dy- 
namics on these sub megaparsec scales. 

However, when MOND is extended to cosmological data sets, 
it seems to comprehensively fail. Apart from the need for dark en- 
ergy, it requires some form of dark matter to fix the dynamics of 
clusters of galaxies and expansion history, as well as the CMB. It 
is therefore obvious that if MOND is to be preserved as a good 
model of galactic dynamics, and at the same time have some mea- 
sure of success in fitting cosmological data, it must do so in coop- 
eration with dark matter. The constraint on this dark matter is that 
it must be hot enough to free stream from galaxies such that it does 
not destroy the good fits to galactic rotation curves. Angus (20(ii) 
showed that several cosmological constraints can be satisfied us- 
ing an 11 eV sterile neutrino that is fully thermalised in the early 
Universe. 

For example, it allows an excellent fit to the angular power 
spectrum of the CMB (reproduced in Fig[T} by virtue of fi^^ for 
an 11 eV sterile neutrino being alm ost identical t o Q.cd m of the 
AC DM model. To complement this. I Angus et alj ^20 id) demon- 
strated that the dark matter density profiles found in MON D of the 
sampl e of groups and clusters of galaxies studied by Ang us et al.l 
( l2008h were consistent with equilibrium configurations of 11 eV 
sterile neutrinos. What this means is that the velocity dispersions 
of the sterile neutrinos required for equilibrium of the density pro- 
files form a phase space density that does not need to exceed the 
Tremaine-Gunn limit. 

Although the phase space densities of sterile neutrinos are 
compatible with static clusters of galaxies, the obvious next step is 
to investigate whether the model can form the observed distribution 
and internal structure of clusters of galaxies in the finite timescale 
allowed by the age of the Universe. Equally as essential, given that 
our model uses hot dark matter, is whether the matter power spec- 
trum is ma tched on the smaller s cales measured by the Lyman alpha 
forest (e.g. lTegmark et al.l2004l) at high redshift. 

Clearly, the only way to test these issu es is w it h cos - 
mological numerical simulations, like those of iNusseij ( |2002|) : 
iKnebe & GibsonI l2004l) : lLlinares et alj f2008l) except that the ini- 
tial conditions will coiTespond, as near as possible, to the cosmo- 
logical model of Angus (2009) - hereafter Angus09. In this pa- 
per we developed a single coarse grid cosmological particle -mesh 
code and then adapted it to solve the modified Poisson equation 
of Milgrom's r ecently proposed QUasi-linear MOND (QUMOND, 
Milgrom 2010). QUMOND is a new incarnation of MOND that is 
more convenient to handle due to the fact that the theory only re- 
quires solutions of a linear differential equation and one non-linear 
algebraic step. It can be derived from an action, thus the conserva- 
tion laws are adhered to. 

Using this QUMOND code, we run a single cosmological sim- 
ulation with as high resolution as possible and compare it with the 
observed distribution of mass on scales of galaxy clusters and be- 
yond. Unfortunately, we do not have the resolution to investigate 
smaller scales that are dominated by baryonic physics. 
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2 QUMOND 

The Aquadratic Lagrangian theory of ISekenstein & MilgromI 

( Il984h produces a modified Poisson equation which must be solved 
numerically for arbitrary geometries. Likewise, QUMOND re- 
quires solution of a modified Poisson equation, but one that is 
slightly easier to implement. Specifically, the ordinary Poisson 
equation for cosmological simulations 



V^$]v =47rG(p-p)/a 



(1) 



is solved to give the Newtonian potential, $jv, at scale factor a, 
from the ordinary matter density p which includes baryons and 
neutrinos. This would also include cold dark matter if there was 
any in our model. The QUMOND potential, <J>, is found from the 
Newtonian potential as follows 



(2) 



where z^(j/) — 0.5 + 0.5^/1 + 4/j/ and y = \/^n /a-oa, with 
Uo being the MOND acceleration constant chosen here to be 
3.6 (kms-^)V"^- 

Essentially, we solve the Poisson equation once, using the or- 
dinary matter density as the source, as would happen in a regular 
N-body code, to find the Newtonian potential (Eq[TJ - incorporat- 
ing periodic boundary conditions. Following this, we take deriva- 
tives of the Newtonian potential and generate a new source density 
(the right hand side of Eq |2j, usually referred to as the ordinary 
matter plus phantom matter density associated with the QUMOND 
field. Then we solve the Poisson equation a second time to find 
the QUMOND potential, "I>. From this, we take the gradient of the 
QUMOND potential to find the gravitational field at each particle's 
location and move our particles accordingly. 

The density assignment follows the standard cloud-in-cell 
technique and from this both Poisson equations (Eqs[T]and|2} are 
solved using Multigrid methods. The densities found from back- 
wards derivation of the potentials gives the source density back 
with impeccable precision at all redshifts and locations in the grid. 
Next comes the important part where the QUMOND source density 
is found. 

Specifically, what happens is the following: assume we want 
the QUMOND source at cell (i,j,k) of a Cartesian grid (x,y,z), then 
we need to define the gravity at various points surrounding it. If we 
use unit length grid cells then 

9x2 ~ 01+1, i,fc ~ 0i,i,fc 

gxi = 4'i,j,k — 4'i-l,j,k 

9y2 = <^i.j+i,fe ~ 
9yi = (t>i,j,k - 

9z2 = <t>i,j,k + l ~ <f>i,j,k 

(?zi = 4'i,j,k — 4>i,j,k-l (3) 

These are the values of the gravitational field at half a cell from 
ii,j,k) in the three orthogonal directions. Similarly, for these six 
points we must find the value of the i/ function. Surrounding the 
point X2 we use 



^1x2 — Yi + l,j,k 



^i,j,k 



(4) 



- <^i + l,i+l,fc + (j^ij + l-.k ~ {4>i + l.j-l,k + (l>i,j-l.k) 
4lJ22 = <^i+l,i,fc + l + 4'i,j,k+l — (0i + l,j,fe-l + 4'i,j,k-l) 



and surrounding xi 

IjJxi ~ 4'i,j,k — <t>i-l,j,k (5) 
4ajyj — (fi,j + l,k + 4>i-l,j+l,k — {4'i,j-l,k + ^i- l.j- 1 ,fe ) 
4aj2j = + 1 + 0!-l,j,fe + l — (0j,j,fe-l + 4'i-l,j,k-l) 

which gives 



(6) 



These are accompanied by Ky2 , , ™d found in a similar 
way. From this we must find 



= '^(^1:2) 

and Vy^ , Vy, , i 



(7) 

Vzi- This finally leaves us with the QUMOND 



source density in cell given by v^^g^^ - v^^^g.:,^ + Vy2gy2 - 

^vi9vi + ^Z9.g z9 ~ t^zi .gzi ■ A good vi sua lisation of the geome try 
can be found in lTiret & CombesI JioOTh or lUlinares et al.lboOSb . 



3 INITIAL CONDITIONS 

The cosmological model we use employs il^^, Q.a, h, 
na)=(0.0443, 0.218, 0.7377, 0.732, 0.955), which means that one 
thermal sterile neutrino is required with a mass near to 11 eV. We 
plot in Fig[T]a comparison between our model's fit to the CMB 
and the ACDM model's. We note that our model requires four ef- 
fective neutrino species (3 active and 1 sterile), and as such the 
*ife fraction expected from big bang nucleosynthesis is 0.261, 
which is roughly 1 — a larger than the measured value (see e.g. 
lManganoet"al]l2010h . The other relic abundances are not nega- 
tively altered and the value for ^Li is sadly not improved. 

We advocate starting simulations before z ~ 100 where 
MOND effects are minimal, and this means producing initial con- 
ditions at these redshifts. As mentioned previously, in the absence 
of a theoretical underpinning of MOND, and given that there is 
no consensus on a relativistic theory, we must produce the initial 
conditions under the ansatz that gravity does not deviate from gen- 
eral relativ ity while z > 100. This allows us to use the COSMICS 
package o f iBertschingeil 11995) without modification. 

The initial power spectrum of the Angus09 model was com- 
puted accurately a t z = 254.1 using the CAMB package of 
iLewis et al. which can be compared (see Fig |2]l to the 

power spectrum of our N-body initial conditions which ar e gen- 
erated using the COSMICS package of lBertschingej ( Il995h . Both 
codes accurately take care of the neutrino distribution function, but 
of course the COSMICS package can only represent the power 
spectrum discretely by distributing particles (in our case 256^) on 
a 257^ grid of length 512 Mpc/h. Each particle weighs 7.8 x 
10" M0. 

The N-body realisation has its amplitude normalised by the 



CMB quadrupole, Qr 



-To 



5C2\2 



, where To = 2.725 K 



is the current temperature of the CMB and C2 is the / = 2 com- 
ponent of the angular power spectrum. This is enforced because 
we cannot use linear theor y in M OND to estimate ag ai z = 0. 
According to Bennett et al., ( l20Ilh . the l-a range of Qrms-ps = 
17 — 45 fiK, so to be highly conservative, we use 18 nK. For com- 
parison, in a AC DAI universe, CJrTns — PS = 18 jj,K corresponds 
to as = 0.68 and the typically used erg = (0.75,0.8,0.9) cor- 
respond to Qrms-PS = (17.5, 18.6, 21.0) p.K respectively. The 
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dependence of the formation of stiucture on this parameter in the 
Angus09 model will be tested in a forthcoming paper. 

Perhaps more significant than this, one can see in Fig|2]that the 
true power spectrum from CAMB is not perfectly reproduced by 
the simulation initial conditions. In particular, the baryonic acous- 
tic peaks and the exact turn off due to the neutrino free stream- 
ing are not replicated exactly. Nevertheless, when the theoreti- 
cally expected power declines due to the neutrino free streaming 
(around A: = 0.1 h/Mpc) the simulated power drops. Another 
obvious feature is that the simulated power at fc > 0.3 h/Mpc 
is somewhat greater than the correct power, however, the power 
is still incredibly small. For instance, for the Angus09 model 
P{k = 0.7 h/Mpc) = 10~^Mpc^//i^, which can be compared 
to a ACDM model (where our £1,^^ is directly swapped for Qcdm) 
which has power at fc = 0.7 h/Mpc of 10~^ Mpc^ /h'^ - some six 
orders of magnitude larger. 



4 THE IMPORTANCE OF MOND FOR STRUCTURE 
GROWTH 

With standard ACDM simulations, we know that from z = 
250 to the perturbations need only grow by a factor of around 
10^ or so at any wavenumber. In contrast, we know that to be vi- 
able models, the sterile neutrino perturbations must increase by 15 
orders of magnitude at large wavenumbers. 

For our test case we took our initial conditions discussed 
in ^ and ran three simulations: one using the code described 
in ij2] and another using only the Newtonian potential i.e. with 
MOND t urned off. We als o ran a simulation using the MLAPM 
code of Knebe & Gibsor j2004l) with the simple fi func ion 



Famaev et alj |2007|) . To be clear, this 



dFamaev & BinnevI |2005| ; 
code does not properly solve the Bekenstein-Milgrom modified 
Poisson equation (it ignores the curl field), but rather only the alge- 
braic relation V<I> = /i(V$) V$jv, where fj, is also an interpolating 
function, similar to the v used in this paper. It is these three simu- 
lations that we discuss in the remainder of the paper. 



4.1 Power Spectrum 

In Fig [3] we plot the z = power spectrum of the particles 
in the three simulations described above along with the measured 
power spectrum of clusters of galaxies from the REFLEX II sur- 
vey (B alaguera- Antolmez et al..,2011.) and also the power spectrum 
of galaxies l^gmark et al.l l2004l) . The final power spectrum of the 
simulation without MOND falls drastically short of the measured 
power spectrum of clusters, and other probes, at all wavenumbers. 
This is the reason why such models with hot dark matter and un- 
tampered gravity are correctly ignored in the literature. 

This power spectrum without MOND can be compared to the 
power spectra found with our QUMOND code and the MLAPM 
simulation where we find that all wavenumbers probed increase in 
power to amplitudes far larger than the power spectrum of galax- 
ies, and, in the case of the QUMOND simulation, far larger than the 
clusters of galaxies. This is potentially problematic, but there exist 
parameters that can be feasibly altered to slow the growth of struc- 
ture - for instance the acceleration constant of MOND decreasing 
with increasing redshift, which is discussed in ^5A\ The serendip- 
itous match between the MLAPM simulation and the power spec- 
trum of clusters may suggest that the gravity missing due to the 
non-inclusion of the curl-field is enough to bring the QUMOND 



simulation into agreement. Alternatively, the disparity may be a re- 
sult of not comparing like for like. For instance, in ACDM, the 
growth of structure is hierarchical and similar on all scales because 
cold dark matter has, by definition, power on all scales. Conversely, 
both MOND and sterile neutrinos have a scale dependence. 

MOND is acceleration dependent and gravity is therefore a 
product of environment, and the sterile neutrinos have a length 
scale below which power is exponentially suppressed. This length 
scale is larger than the diameter of galaxies and thus sterile neutri- 
nos have merely a shepherding influence on the baryons that build 
galaxies. For example, the sterile neutrinos collapse (together with 
the baryons) into filamentary structures, with clusters of galaxies 
forming at the nodes of these filaments. At this point, galaxies can 
only form if a perturbation from tidal torques or some other shock 
to the primordial gas can create a perturbation which can grow un- 
der MOND. Obviously this cannot occur efficiently at the nodes 
where the galaxy clusters form because the external gravitational 
influence would diminish the effect of MOND and galaxies cannot 
form without boosted gravity (whether be it from cold dark matter 
or modified gravity). 

For these reasons, an accurate comparison with the galaxy 
power spectrum will have to wait until the software of MOND cos- 
mological simulations improves dramatically. 

To complement the comparison between structure formation 
with MOND and without, we also show in Fig|4]a two dimensional 
projection of the density contrast in the two different scenarios. 
Specifically plotted are the particles in a 50 Mpc/h slice along 
the line of sight, where each particle is given a colour to indicate 
the local three-dimensional density. The colour contrasts in the two 
figures are different, the maximum contrast with MOND at this res- 
olution is many thousand times the cosmic mean, without MOND it 
is barely 50% . Therefore, any successful cosmology based on ster- 
ile neutrinos with an 11 eV mass needs a modified gravity theory 
which boosts perturbations, especially on small scales. 



5 GALAXY CLUSTERS 
5.1 Mass function 

Despite a comparison with the galaxy power spectrum be- 
ing unfeasible due to the finite resolution (each particle weighs 
7.8 X IO'^^Mq), it is important that the structures formed in the 
simulation at 2: = resemble the observed distribution and den- 
sit ies of clusters of gala xies - which can be adequately resolved. 
In lPierpaoUetalJ ( I2OOII) , hereafter PSWOl, their Fig 6 shows the 
observed cumulative number distribution of X-ray clusters above a 
specific temperature. The expected number of X-ray clusters with 
a temperature larger than 4.5 ± 0.5 keV and 8.5 ± 0.5 kcV is 
107 - 200 and 14 - 28 respectively in a (512 Mpc/hf volume. 
Upon the z = Q output o f our simulations, we ran the halo finder of 
iKnebe & Gibsorj |2004 ) and discarded all halos with less that 1000 
particles, so that we could adequately resolve the density profiles 
of all 155 cluster mass halos. We show in Fig |5(a)| the cumulative 
number of halos in our (512 Mpc/h)^ volume as a function of 
mass enclosed within 1 Mpc (0.732 Mpc/h). All these 155 halos 
had masses within 1 Mpc of greater than lO^'*i\f0. 

In these simulations, we are frequently required to compare 
the mass distributions of our MONDian halos, with those that a 
Newtonist would derive. Obviously our halos are composed of 
baryons and sterile neutrinos and have a specific gravity profile de- 
pending on the mass profile - that is, the MONDian mass profile. 
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Similarly, this gravity profile observed by a Newtonist would cor- 
respond to another, Newtonian mass profile. In QUMOND, there 
is a simple algebraic relationship between a spherical MONDian 
mass profile in MONDian gravity and the Newtonian gravity of that 
MONDian mass profile, and therefore, the Newtonian mass profile. 
The Newtonian gravity of the MONDian mass profile Adrn{r) is 
naturally (;jv(r) — GMm{r)/r^ and the MONDian gravity profile 
of the MONDian mass distribution is gM(r) — h'[gN{r))gN{r). 
The final step is to say that this MONDian gravity profile corre- 
sponds directly to a (fictitious) Newtonian mass profile i.e. the mass 
profile a Newtonist would derive from a dynamical probe of some 
kind, and this is simply Afjv(r) = gAi{r)r^ /G. Obviously this 
means the radius at which you compare your mass estimates in the 
two gravity theories can make a huge difference to your mass ratio, 
which is worth keeping in mind. 

Therefore, at the radius of 1 Mpc which we chose to 
compare with some results from the literature, a MOND mass 
of IO^^Mq corresponds to a Newtonian gravity of (7jv(r = 
IMpc) = 4.42 X 10"^ X 10^''/(10'')^ = 0.442( kms-i)Vpc. 
Now since the MOND acceleration constant in these units is 
Uo = 3.6 (kms"^)Vpc we find that x = 0.442/3.6 = 0.123, 
i^iy) = 3.4 and therefore the Newtonian mass is 3.4 x lO^M©. 
Thi s large a mass is typica lly found in clusters like A133 and A383 
(see lVikhlinin et al.l2006h at 4.1 to 4.8 keV respectively. Similarly, 
we expect a MOND mass of 4 x lO^^M© in clusters like the 8-9keV 
A 2029 and A 2390. For these typically observed MOND masses, 
we use the range of expected clusters from PSWOl as a gauge of 
whether our simulation forms the correct number of halos. 

It would appear from Fig |5(a)| that we form roughly the correct 
order of magnitude of clusters hotter than 4.5 keV, but form far too 
many high mass clusters. This contradicts the naive expectation that 
galaxy clusters cannot be formed with MOND and 11 eV sterile 
neutrinos. It is in fact reassuring that there is a spectrum of cluster 
masses formed and that more low mass clusters are formed at the 
expense of larger clusters and that all the matter does not merely 
end up in a single cell at the centre of the simulated grid. 

We have included a similar test in Fig |5(b)| in which is plotted 
the number density of clusters per (Mp c/hY and per int erval of 
logioA'hoo- The data points come from^ Rines et al.l 12008 ). here- 
after RDN08, where the filled circles are measurements using a 
scaling relation between X-ray luminosity and the virial mass and 
the open diamonds are usin g the caustic method o f iDiaferidd 19991) 
(see also ISerra et al.ll201 ih to estimate M2oo- Our estimation us- 
ing the halos found in our QUMOND simulation is computed in 
a similar way to the aforementioned comparison with the data of 
PSWOl. Again we begin by initially finding the MONDian mass 
profile of all the clusters and then convert it to the Newtonian mass 
which would provide the same gravity as a function of radius. From 
the Newtonian mass we derived the Newtonian A/200 which is the 
enclosed mass at the radius where the density falls to 200 times 
the critical density. Given that a Newtonian virial mass of 10^^ Mq 
at several Mpc would correspond roughly to a MONDian mass of 
lO^*Af0, our resolution precludes us from comparing with all but 
that last data point shown in RDN08. 

Our model appears to demonstrate the same trend with the 
RDN08 data as it does with the PSWOl data, except that the 
RDN08 data probes somewhat smaller cluster masses than the 
PSWOl sample. By probing cluster temperatures up to 9 keV, 
the PSWOl sample is probing Newtonian virial masses up to 
2 X IO^^AIq {logioA'l2oa ~ 15.3), which is marginally larger 
than the logioM2oo = 15.1 limit of RDN08, which by contrast 
probes far lower virial masses than the PSWOl survey does. This 



is the reason why our line goes through the last two data points in 
Fig |5(b)| and yet over estimates the number of very massive clusters 
(T = 8.5 keV) in Fig [5(a)l Our hypothesis is that if the RDN08 
study was extended to io(;io A/200 = 15.3, our current model would 
overestimate the numbers per mass bin. 

Comparing Figs |5(a)| and |5(b)| it is apparent that we have 
agreement with the moderately massive clusters (T — 4.5 keV) in 
Fig |5(a)| only because we produce many times more T > 8.5 keV 
clusters than desired. Therefore, Fig |5(b)| demonstrates that we 
severely underpredict the number of T = 4.5 keV clusters. This 
currently is a resolution issue that will be addressed by future sim- 
ulations. A point worth noting is that the mass function per mass 
bin as plotted in Fig |5(b)| is flat for all halo masses. However, if so 
much mass was not locked up in these extremely massive MOND 
halos with > 5 x IO^^A/q within 1 Mpc (in fact, there is at least 
5 X 10^^ Mq locked up in these), then there could be 500 more 
10^'* A/0 halos produced in principle. 

The task now is tuning the available free parameters in order 
to fit the measured mass function. These include the normalisation 
of the initial power spectrum, the i^-function as well as two con- 
siderably more interesting prospects. The first is how the MOND 
acceleration constant varies with redshift and the second is how the 
background expansion of a QUMOND cosmology might vary from 
the FRW family of models. The first case is almost unconstrained 
since no detailed studies of the dynamics of galaxies exist at any 
significant redshift. 

Here we make only the hypothesis that the increased energy 
density of dark energy with increasing scale-factor is linked (per- 
haps directly) to the increased prevalence of QUMOND fields as 
the Universe becomes more sparse. If accelerations are much larger 
than Qo a.1 z > 100 then there would be no QUMOND fields, but 
as the Universe expands and accelerations drop below ao (or what- 
ever value Uo takes at a given redshift) then these QUMOND fields 
would begin to increase. This is qualitatively the same trend as dark 
energy takes and is a suggestion that will be expanded upon in a 
forthcoming paper where we will constrain models of the expan- 
sion history of the Universe with supemovae and gamma-ray burst 
data. 

In fact, using a full Bayesian appraoch. iDiaferio et al.l Jio 1 ih 
have shown that data from supernovae and gamma-ray bursts can 
be described by the KCDM model as well as other more exotic 
expansion histories. For now it is enough to remark that a different 
expansion history could be key to the issues mentioned above, with 
reference to Figs |5(a)| and |5(b)| for the following reasons: if the for- 
mation of the over-abundance of very massive clusters (> 10^^ A/© 
in terms of their MOND mass) that we see in our simulations is due 
to the merging of many subhalos or alternatively due to monolithic 
collapse, then more rapid expansion at a certain key stage could 
halt the formation of these superclusters and sustain the abundance 
of lower mass halos. 

Most likely, the expansion history of the Universe, a varying 
ao, the i^-function and the normalization of the power spectrum 
must be traded off against each other. In addition, our resolution 
limit of roughly 10^^A/q within 1 Mpc must be greatly extended 
in order to solidify these claims by checking the formation of clus- 
ters down to small grou ps of galaxies, whic h have been shown in 
lAngus et all 1 I2OO8I) and lAngus et all ( I2OI0I) to host dense sterile 
neutrino halos - IQ^^Mq within 0.2 Mpc - and prove to be the 
strictest test of the Angus09 model. 

Although it would appear now to be within reasonable doubt 
that we could fit the distribution - at least the high temperature end - 
of massive clusters of galaxies, this in no way gives any information 
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about the formation and distribution of galaxies. As we discussed 
in ij4] the formation process of galaxies is entirely dependent on 
perturbations to the distribution of cold gas and the accretion of 
more gas aided by MONDian gravity. The sterile neutrinos will 
have no influence except presumably as distant tidal torques. 

5.2 Density profiles 

Another thing that is uncertain about galaxy clusters formed 
from sterile neutrinos under M ONDian gravity is their internal 
structure. In lAngus et al.l ( 12010. ) we extracted the density profiles 
of the sterile neutrinos required to fit the mass profiles of X-ray 
clusters measured from the density and temperature of the intra- 
cluster medium. After merging the sterile neutrino densities with 
the small contribution from the X-ray gas, we show in Fig |6(a)| 
that the two largest clusters (A 2029 and A 2390) both had log- 
arithmic density slopes over the range 0.1 to 1.0 Mpc of around 
^1"^ ~ —2.0, and certainly greater than -2.5. We plot alongside 
this the density profiles for four typical halos from our simulation 
with similar masses. Our clusters have outer density slopes that are 
significantly shallower than —2. The shallower densities are due 
to a lack of resolution and to avoid the conclusion that this is a 
generic feature of sterile neutrino halos emerging from cosmolog- 
ical simulations, we also plot in Fig |6(b)| the density profile of our 
most massive halo out to 25 Mpc, where resolution is certainly not 
a problem. We have overplotted logarithmic slopes of -2 and -1.8, 
and the latter seems to resemble the outer density profile over an 
order of magnitude in radius. 

In addition to this, baryonic adiabatic contraction in clusters 
has been shown by Gnedin et al. (2004) to increase the steepness 
of the density profiles of cluster dark matter halos in a A.CDM 
universe, so there remain several reasons to be optimistic about the 
densities of cluster size halos at intermediate radii. Were we to have 
perfect resolution, we would expect our halos to resemble the best 
fits to cluster of galaxy dark matter halos in MOND, at least down 
to 100 kpc. That is unless there is a scale dependence on cluster ha- 
los such that the most massive clusters have significantly shallower, 
or steeper, outer density profiles. 

Related to the discussion in i|5.1| it may be the case that the 
steepness of the halos is dependent on the redshift of formation 
which would be different if the MOND acceleration constant was 
lower at higher redshift. What is certain is that the formation red- 
shift of all clusters will drop if the MOND acceleration constant is 
lower at high redshift, but the effect on the steepness of the density 
profiles is less clear. In CDM simulations for example, the concen- 
tration parameters of halos are observed to decrease with increasing 
mass (hence decreasing formation redshift) i. e. smaller halos form 
earlier and are therefore more dense (see e.g. iBullock et alll200lh . 
Presumably this result holds also in MOND, so delaying the for- 
mation to lower redshift would decrease the concentration of the 
sterile neutrino halos. Increasing the amplitude of the initial fluc- 
tuations (through the CMB quadmpole) could act in the opposite 
direction, so perhaps tuning between the two can provide the cor- 
rect balance between number density of halos and steepness of the 
density profiles of individual clusters. 

Equally importan t, the inner profiles of the clusters studied 
m lAngus et alj jlQl& i were found to reach densities in excess of 

Mqpc~^ at 10 kpc. This was the case for all but the sub 
1 keV groups of galaxies. At the moment, our simulations are un- 
reliable below 1000 kpc/h and cannot probe such small radii, but 
it will be vital in future studies to check whether the densities of 
the sterile neutrino halos can continue to rise all the way to 10 kpc 



and reach the relatively high densities found there in the real cluster 
halos. 

The reason we cannot simply run a simulation with a 
16 Mpc/h box and immediately find the central densities of the 
clusters is because that small a box cannot incorporate the influ- 
ence of the large scale structure on the local environment, which is 
of crucial importance in MOND cosmological (and ordinary) simu- 
lations. Ideally, one should not use a box smaller than we have used 
here and so a refinement strategy and more than 256'' particles will 
be necessary. In princip le, a resimulation technique, like that used 
bv lSpringel et al.l J200l[) , could be invaluable. 



6 PAIRWISE VELOCITIES 

The r elative veloc ity of the two clusters that comprise the bul- 
let cluster telowe et a l. 2004, 2006; Bradac et al. 2006) is a topic 
of much interest basically because, at first order, it appears to 
be larger than one could naturall y expect in a ACDM univers e 
(see e.g. [M arkevitch et al. 20021; iMarkevitch & VikhlininI l2007h . 
iMastropietro & Burkerli 02008) have recently shown via hydrody- 
namic simulations that the infall velocity of the clusters must be 
3000 kms^^ within 2-3i?200 in order to reprod uce the observed 
features of the various shocks. With this in mind, iLee & Komatsul 
(f2010) have simulated a large cosmological volume (27 Gpc^ /h^) 
to infer the statistical probability of bullet cluster like events, which 
they showed to be incredibly low (between 3.3 x 10^^^ and 
3.6 X 10"^). 

To estimate the likelihood of creating bullet cluster like events 
in the Angus09 cosmological model, we show the z = (not 
z = 0.3) relative velocities between all pairs of halos within a 
sphere of 10 Mpc from each other, which is in general between 
3-6^200. We find that out of the 155 x (155 - 1) /2 ^ 12000 po- 
tential pairs, there are less than 100 pairs within 10 Mpc of each 
other and out of these we find 5 with relative velocities exceed- 
ing 3000 kms~^ (see Fig|7j. These velocities will increase with 
decreasing separation. 

This number is inclusive of all bullet cluster sized halos since 
the sub cluster was observed to have a temperature of 4.5 keV, 
which is roughly the limit of our resolution. Although the relative 
velocities will b e lower at z = 0.3 it still se ems plausible that, just 
as was found bv lAngus & McGaughl ( l2008h . MOND has the ability 
to produce very swift collisions between clusters. With better reso- 
lution it might be interesting to check the contrary problem which 
is if MOND would produce too many bullet cluster like systems. 



7 CONCLUSION 

Here we introduced a new Particle-Mesh code that accurately 
solves the modified Poisson equation of Milgrom's QUMOND the- 
ory. We presented a set of initial conditions for the Angus09 cos- 
mological model and computed the power spectrum to demonstrate 
how modified gravity (in particular MOND) is absolutely essential 
to create large scale structure in a hot dark matter universe. 

We ran a simulation from z — 250 to 2 = and identified 
that we form roughly the correct order of magnitude number of 
massive clusters of galaxies hotter than 4.5 keV and form a rea- 
sonable spectrum of masses. Having said that, we considerably un- 
derpredict the number of T = 4.5 keV clusters and overpredict 
the number of T = 8.5 kcV clusters. We suggest that the former 
problem is a resolution issue that will be dealt with by forthcoming 
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Figure 4. Shows the distribution of particles in a (512 Mpc/h) X (512 Mpc/h) X (50 Mpc/h) slice, where the particles are colour coded according to 
their local 3D density. The left hand panel is a simulation without MOND and the right hand panel is with MOND. The contour levels are very different in 
the two figures: the left hand panel's extreme density contrast is only 50% greater than the cosmic mean, whereas the right hand panel's is many thousands of 
times the cosmic mean. 
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Figure 5. (a) In the left hand panel we show the cumulative number of halos in a (512 Mpc/h)^ volume as a function of mass enclosed within 1 Mpc. 
Against this curve we have plotted the number of halos hotter than 4.5 keV and 8.5 keV in such a volume as measured by Pierpaoli, Scott and White (2001). 
(b) In the right hand panel we plot the number density of halos per (Mpc/h)'^ and per logarithmic interval of mass plotted against logio of A/200- The data 
points come from Rines, Diaferio & Natarayan (2008). The filled circles are measurements using a scaling relation between X-ray luminosity and the virial 
mass, whereas the open diamonds are using the caustic method of Diaferio (1999) to estimate A/200- The line is using the halos found in our QUMOND 
simulation. 



simulations, but that some variation of the available free parame- 
ters will be necessary to match the number of very rich clusters of 
galaxies. The two key variables are how the MOND acceleration 
constant decreases with increasing redshift and how the expansion 
history of our model differs from the FRW prediction. 

Our resolution was insufficient to compare the density profiles 
of our halos with those required to fit X-ray clusters in MOND 
jAngus et alj 12010.) . However, we showed that the logarithmic 
slope of the largest clusters in our simulation, which extended to 
beyond 25 Mpc and hence well into the resolved range, was com- 
patible with the observed halos. 

To follow up the work of I Angus & McGaugl] ( 1200 Sb . which 
suggested that the large relative velocity of the bullet cluster was 
a natural result of MOND, we found the numbers of pairs of clus- 
ters within 10 AIpc of each other and used their relative velocitie s 
to compare with the conclusion of iMastropietro & Burkerj fcOOSl) . 
They showed that the relative velocity of a pair needed to be above 
a certain threshold velocity to generate the observed bow shock and 



other features of the bullet cluster system. We found that large rela- 
tive velocities are indeed a natural product of a MOND cosmology. 

In the future, it will be vital to amend the code to resolve struc- 
tures down to below 100 kpc and even 10 kpc. Only then will it be 
clarified whether or not the cluster halos formed in cosmological 
N-body simulations with MOND are fully compatible with the ob- 
served mass discrepancy in MONDian clusters. 
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Figure 6. (a) In the left hand panel the lines plotted are the density profiles of cluster sized halos with a similar mass to the rich clusters A 2029 and A 2390, 
whose MONDian density profiles are also plotted as the triangle and filled circle symbols. We note that the radiis probed here are below the resolution limit of 
the simulation, (b) In the right hand panel we plot the density profile of the largest cluster we form which extends well into the resolved range and overplotted 
are logarithmic density slopes of -2 and -1.8 with dashed linetype. 
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Figure 1. Shows the CMB angular power spectrum for our cosmological 
model (blue line), compared with the KCDM model (red line). The data 
points come from WMAP year 7 (black), ACT (turquoise) and ACBAR 
(green). 
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Figure 2. Shows the power spectrum of our N-body realisation (dashed) at 
2 = 254.1, the expected power spectrum from the CAMB package (solid) 
for that model as well as the power spectrum for the AC DM model (dotted 
line). 
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Figure 3. Shows the measured power spectrum of galaxy clusters at 2 = 
from the REFLEX II survey (filled circles) and for reference the galaxy 
power spectrum (open triangles) from Tegmark et al. (2004). We also plot 
the power spectra of particles from our three simulations: QUMOND us- 
ing our code (solid), algebraic MOND using MLAPM (dotted) and without 
MONO (dashed). 

Martin Feix, Benoit Famaey, Gianfranco Gentile and Kurt van der 
Heyden. 
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